Przyglądnijmy się teraz innej klasie problemów obliczeniowych, tak zwanemu problemowi Stokes'a. Równania te opisują przepływ cieczy [1], [2].
Spróbujmy najpierw zrozumieć sens tych rówńań. W jaki sposób można opisać przepływ cieczy na obszarze dwuwymiarowym? Załóżmy że szukamy rozwiązania problemu Stokesu (problemu przepływu cieczy) na obszarze w kształcie kwadratu \( \Omega = [0,1]^2 \). W każdym punkcie \( (x,y)\in \Omega \) naszego obszaru chcemy znać prędkość przepływu cieczy, która jest wektorem i każdemu punktowi przypisuje dwie współrzędne pola prędkości (prędkości w kierunku dwóch osi układu współrzędnych w tym punkcie) \( \Omega \ni (x,y) \rightarrow {\bf u(x,y) }=(u_1(x,y),u_2(x,y))\in R^2 \) (gdzie pogrubione \( {\bf u } \). oznacza wektor dwóch składowych) oraz pole skalarne ciśnienia, które każdemu punktowi przypisuje jedną wartość ciśnienia w tym punkcie \( \Omega \ni (x,y) \rightarrow p(x,y)\in R \).
Mamy więc trzy niewiadome w każdym punkcie, dwie składowe prędkości oraz jedną składową ciśnienia. Potrzebujemy więc co najmniej trzech równań żeby rozwiązać nasz problem. Równanie Stokesa w ogólnej formie zapisuje się w postaci
\( - \Delta {\bf u } + \nabla p = f \).
Pamiętając że pole prędkości ma dwie składowe, oraz przypominając sobie co oznacza Laplasjan oraz gradient, dostajemy dwa równania w tak zwanej formie silnej
\( -\frac{\partial^2 u_1(x,y)}{\partial x^2} -\frac{\partial^2 u_1(x,y)}{\partial y^2} + \frac{\partial p(x,y)}{\partial x } = f_1(x,y) \)
\( -\frac{\partial^2 u_2(x,y)}{\partial x^2} -\frac{\partial^2 u_2(x,y)}{\partial y^2} + \frac{\partial p(x,y)}{\partial y } = f_2(x,y) \)
Dodatkowe trzecie równanie dostajemy zakładając pewne dodatkowe własności cieczy której przepływ chcemy obliczyć. Tradycyjnie, w równaniu Stokesa zakłada się że ciecz jest nieściśliwa, oraz że nie mamy źródeł ani ubytków cieczy w obszarze. Matematycznie zapisuje się ten warunek w postaci tak zwanej dywergencji, piszemy konkretnie że dywergencja jest zerowa \( \displaystyle{div} {\bf u}=0 \) gdzie \( \displaystyle{div} {\bf u} = \nabla \cdot {\bf u} = \left( \frac{\partial }{\partial x}, \frac{\partial }{\partial y } \right) \cdot \left(u_1,u_2 \right) = \frac{\partial u_1(x,y)}{\partial x} +\frac{\partial u_2(x,y)}{\partial y } =0 \)
Mamy więc następujący układ równań do rozwiązania
\( -\frac{\partial^2 u_1}{\partial x^2} -\frac{\partial^2 u_1}{\partial y^2} + \frac{\partial p}{\partial x } = f_1 \)
\( -\frac{\partial^2 u_2}{\partial x^2} -\frac{\partial^2 u_2}{\partial y^2} + \frac{\partial p}{\partial y } = f_2 \)
\( \frac{\partial u_1}{\partial x} +\frac{\partial u_2}{\partial y } =0 \)
Żeby rozwiązać równanie różniczkowe cząstkowe które spełnione jest na danym obszarze, musimy dodatkowo określić co dzieje się na brzegu tego obszaru. W ten sposób sprecyzujemy problem który chcemy rozwiązać. Każde inne warunki brzegowe dają nam inny problem i w konsekwencji inne rozwiązanie. W naszym przypadku wprowadzamy tak zwane warunki brzegowe Dirichleta na prędkość przepływu, które to mówią nam jakie są wartości pola prędkości na brzegu obszaru.
Załóżmy że nasz obszar jest to zatoka na brzegu rzeki, która płynie jednostajnie z lewej strony na prawą, wzdłuż górnej granicy naszego obszaru \( \Omega \). Zakładamy że prędkość przepływu na górnym brzegu obszaru wynosi
\( {\bf u}(x,y) = (1,0) \textrm{ dla } x \in (0,1), y=1 \).
Na pozostałej częśi brzegu zakładamy że prędkość przepływu wynosi zero
\( {\bf u}(x,y) = (0,0) \textrm{ dla } x \in \{0,1\}, y\in (0,1) \cup x\in(0,1),y=0 \).
Mamy więc określone warunki brzegowe definiujące wartości pola prędkości na brzegu obszaru. Warunki takie nazywamy warunkami brzegowymi Dirichleta. W problemie Stokesa obliczamy jeszcze ciśnienie, musimy więc przyjąć jakieś założenie o ciśnieniu. Przyjmujemy warunek iż sumarycznie ciśnienie na całym obszarze jest równoważone (ciśnienie próbkowane w różnych miejscach poprzez przemnożenie przez funkcję testującą oraz odcałkowanie, daje w sumie zero)
\( \int_{\Omega} p(x,y) dxdy = 0 \).
Możemy teraz skonstruować tak zwane sformułowanie słabe naszego problemu, które będzie nadawać się do symulacji metodą elementów skończonych.
Sformułowanie słabe uzyskuje się w sposób następujący. Całkujemy i przemnażamy nasze równanie przez wybrane funkcje zwane funkcjami testującymi których używać będziemy do uśredniania naszego równania w obszarze na którym funkcje te są określone. Główna różnica w stosunku do problemu projekcji jest taka że mamy teraz układ trzech równań, potrzebujemy więc trzech funkcji testujących, po jednej do każdego równania, oznaczamy je \( v_1(x,y), v_2(x,y), q(x,y) \). Przemnażamy równania i całkujemy po obszarze
\( -\int_{\Omega} \frac{\partial^2 u_1(x,y)}{\partial x^2}v_1(x,y) dxdy -\int_{\Omega} \frac{\partial^2 u_1(x,y)}{\partial y^2}v_1(x,y) dxdy + \int_{\Omega} \frac{\partial p(x,y)}{\partial x }v_1(x,y) dxdy = \int_{\Omega} f_1(x,y) v_1(x,y) dxdy \)
\( -\int_{\Omega} \frac{\partial^2 u_2(x,y)}{\partial x^2}v_2(x,y) dxdy -\int_{\Omega} \frac{\partial^2 u_2(x,y)}{\partial y^2}v_2(x,y) dxdy + \int_{\Omega} \frac{\partial p(x,y)}{\partial y } v_2(x,y) dxdy= \int_{\Omega} f_2(x,y) v_2(x,y) dxdy \)
\( \int_{\Omega} \frac{\partial u_1(x,y)}{\partial x}q(x,y) +\int_{\Omega} \frac{\partial u_2(x,y)}{\partial y }q(x,y) dxdy =0 \)
W podobny sposób jak miało to miejsce w przypadku problemu projekcji bitmapy, każdy wybór funkcji testującej \( v(x,y) \) służącej do uśredniania naszego problemu w obszarze w którym funkcja testująca jest określona, daje nam jedno równanie. Różne wybory funkcji testującej \( v(x,y) : v_1(x,y),v_2(x,y),q(x,y) \) dadzą nam więc różne równania
\( \int_{\Omega} \frac{\partial u_1(x,y)}{\partial x}\frac{\partial v_1(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial u_1(x,y)}{\partial y}\frac{\partial v_1(x,y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial u_1(x,y)}{\partial n} v_1(x,y) dS + \int_{\Omega} \frac{\partial p(x,y)}{\partial x }v_1(x,y) dxdy = \int_{\Omega} f_1(x,y) v_1(x,y) dxdy \)
\( \int_{\Omega} \frac{\partial u_2(x,y)}{\partial x}\frac{\partial v_2(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial u_2(x,y)}{\partial y}\frac{\partial v_2(x,y)}{\partial y} dxdy + \int_{\partial \Omega} \frac{\partial u_2(x,y)}{\partial n} v_2(x,y) dS + \int_{\Omega} \frac{\partial p(x,y)}{\partial y }v_1(x,y) dxdy = \int_{\Omega} f_2(x,y) v_2(x,y) dxdy \)
\( \int_{\Omega} \frac{\partial u_1(x,y)}{\partial x}q(x,y) +\int_{\Omega} \frac{\partial u_2(x,y)}{\partial y }q(x,y) dxdy =0 \)
gdzie po całkowaniu przez części pojawiają się nowe człony - całki brzegowe liczone w kierunku prostopadłym do brzegu.
W naszym problemie Stokesa związanym z obliczaniem przepływu w zatoce sąsiadującej z rzeką zakłada się że siły \( f_1(x,y)=f_2(x,y)=0 \). W warunku brzegowym mamy założone że prędkość na górnej krawędzi wynosi
\( {\bf u}(x,y) = (1,0) \textrm{ dla } x \in (0,1), y=1 \). Innymi słowy wprowadzamy tutaj warunek brzegowy Dirichleta na prędkość przepływu.
Ze sformułowaniem tym wiąże się jeden techniczny problem. Otóż pole przepływu na górnym brzegu jest niezerowe (skierowane w prawo) natomiast na bocznych brzegach (na lewym i prawym brzegu) pole przepływu przyjmujemy równe zero. Mamy więc przeskok z zera do jeden na pierwszej współrzędnej pola prędkości. Przeskok taki po pierwsze nie jest fizyczny (nie może być tak że woda w jednym punkcie nie płynie a zaraz obok płynie całkiem szybko, przejście z "nie płynięcia" do "płynięcia" odbywa się w gładki regularny sposób. Po drugie, jeśli chcielibyśmy rozwiązanie numeryczne przybliżać za pomocą metody elementów skończonych, rozwiązanie takie jest przedstawiane jako kombinacja liniowa wielomianów ciągłych. Nie ma takich wielomianów które przeskakują z zera do jeden. Postawienie problemu w taki sposób spowoduje powstanie osobliwości numerycznych w lewym górnym i prawym górnym brzegu obszaru. Żeby tego uniknąć, w tym przypadku wykonuje się następujący trick. Przyjmuje się
\( g(x,y) = 1 \textrm{ dla } x \in (\epsilon, 1-\epsilon), y=1 \)
\( g(x,y) = \frac{x}{\epsilon} \textrm{ dla } x \in (0,\epsilon), y=1 \)
\( g(x,y) = \frac{(1-x)}{\epsilon} \textrm{ dla } x \in (1-\epsilon, 1), y=1 \)
\( {\bf u}(x,y)=(g(x,y),0) \textrm{ dla } x \in (0,1), y=1 \)
Wówczas faktycznie w rogach obszaru mamy zerową prędkość i nie ma osobliwości numerycznych w obliczeniach metodą elementów skończonych.
Dokonujemy następnie tak zwanego "przesunięcia" naszego problemu. Poszukiwane rozwiązanie dzielimy na dwie składowe
\( {\bf u}(x,y) = (h(x,y)+w(x,y),u_2(x,y)) \)
gdzie \( h(x,y)=g(x,y)(1-y) \) jest tak zwanym rozszerzeniem naszego niezerowego warunku brzegowego Dirichleta na cały obszar \( \Omega \).
Naszego przesunięcia dokonujemy więc tak naprawdę jedynie względem pierwszej składowej pola prędkości (ponieważ druga składowa jest równa zero). Wówczas pierwsze równanie
zmienia się
\( \int_{\Omega} \frac{\partial (h(x,y)+w(x,y))}{\partial x}\frac{\partial v_1(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial (h(x,y)+w(x,y))}{\partial y}\frac{\partial v_1(x,y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial (h(x,y)+w(x,y))}{\partial n} v_1(x,y) dS + \int_{\Omega} \frac{\partial p(x,y)}{\partial x }v_1(x,y) dxdy = 0 \)
Człon z zadaną funkcją \( h(x,y) \) jest znany i możemy przenieść go na prawą stronę
\( \int_{\Omega} \frac{\partial w(x,y)}{\partial x}\frac{\partial v_1(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial w(x,y)}{\partial y}\frac{\partial v_1(x,y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial w(x,y)}{\partial n} v_1(x,y) dS + \int_{\Omega} \frac{\partial p(x,y)}{\partial x }v_1(x,y) dxdy = \\ = \int_{\Omega} \frac{\partial h(x,y)}{\partial x}\frac{\partial v_1(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial h(x,y)}{\partial y}\frac{\partial v_1(x,y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial h(x,y)}{\partial n } v_1(x,y) dS \)
Wszystkie człony znajdujące się po prawej stronie potrafimy obliczyć. Z koleji niewiadomą \( w(x,y) \) wprowadzoną na lewej stronie możemy z powrotem oznaczyć jako \( u_1(x,y) \) musimy jednak pamiętać po rozwiązaniu naszego problemu Stokesa że należy do niej później dodać z powrotem człon z rozszerzeniem warunku brzegowego.
Wprowadzamy oznaczenia na naszą lewą i prawą stronę
\( a({\bf u}, {\bf v}) =\int_{\Omega} \nabla {\bf u} : \nabla {\bf v}dxdy \)
gdzie \( : \) oznacza sumowanie wszystkich składowych
\( =\int_{\Omega} \frac{\partial u_1}{\partial x}\frac{\partial v_1}{\partial x}dxdy + \int_{\Omega} \frac{\partial u_1}{\partial y}\frac{\partial v_1}{\partial y}dxdy + \int_{\Omega} \frac{\partial u_2}{\partial x}\frac{\partial v_2}{\partial x}dxdy + \int_{\Omega} \frac{\partial u_2}{\partial y }\frac{\partial v_2}{\partial y }dxdy \)
\( b(p, {\bf v}) = \int_{\Omega} p div {\bf v} dxdy \)
\( = \int_{\Omega}p \frac{\partial v_1}{\partial x} dxdy+ \int_{\Omega} p \frac{\partial v_2}{\partial y } dxdy \)
\( f({\bf v} )= \int_{\Omega} \frac{\partial h(x,y)}{\partial x}\frac{\partial v_1(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial h(x,y)}{\partial y}\frac{\partial v_2(x,y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial h(x,y)}{\partial n } v_1(x,y) dS \)
Nasz problem można teraz zapisać w następującej postaci
\( a({\bf u},{\bf v}) + b(p,{\bf v }) = f({\bf v}) \quad \forall {\bf v} \)
\( b(q,{\bf u } ) = 0 \quad \forall q \)
Zauważmy że zastąpiliśmy układ trzech równań układem dwóch równań. Nasze dwa równania są "zakodowane" w pierwszym równaniu. Teraz, możemy przyjąć funkcje testujące \( {\bf v }=(v_1,0) \) żeby dostać pierwsze równanie, oraz \( {\bf v }=(0,v_2) \) żeby dostać drugie równanie.
Mając w ten sposób opisany nasz problem, możemy teraz przyjąć że nasze poszukiwane rozwiązanie \( ({\bf u },p)= (u_1(x,y),u_2(x,y),p(x,y)) \) ma człony prędkości równe zero na brzegu \( u_1(x,y)=u_2(x,y)=0\textrm{ dla }(x,y)\textrm{ na }\partial \Omega \). Podobnie możemy przyjąć że nasze funkcje testujące, są również równe
zero na brzegu \( v_1(x,y)=v_2(x,y)=0\textrm{ dla }(x,y)\textrm{ na }\partial \Omega \)
Określamy teraz funkcje używane do aproksymacji rozwiązania oraz do testowania.
Precyzujemy jak wyglądać będzie nasza grupa elementów, oraz nasze funkcje bazowe B-spline rozpięte na elementach, podając wektor węzłów wzdluż osi \( x \) oraz wektor węzłów wzdłuż osi \( y \). Odsyłamy tutaj do rozdziału trzeciego w którym opisany jest sposób definiowania funkcji bazowych za pomocą wektorów węzłów i wielomianów B-spline.
Wzdłuż osi \( x \) wprowadzamy wektor węzłów [0 0 0 1 2 3 4 4 4], podobnie wzdłuż osi \( y \) wprowadzamy wektor węzłów [0 0 0 1 2 3 4 4 4].
Uzyskaliśmy dwie bazy jednowymiarowych funkcji B-spline \( \{ B_{i,2}(x) \}_{i=1,...,6 } \) oraz \( \{B_{j,2}(y)\}_{j=1,...,6 } \). Następnie utworzymy z nich iloczyny tensorowe \( B_{i,j;2}(x,y)=B_{i,2}(x)B_{j,2}(y),i,j=1,...,6 \).
Zauważmy że nasz obszar \( \Omega \) rozpięty jest na kwadracie \( [0,1]^2 \), podobnie jak nasze 6*6=36 funkcji bazowych co wynika z definicji wektora węzłów [0 0 0 1 2 2 3 4 4 4].
Przyjmujemy teraz naszą aproksymację pól prędkości
\( u_1(x,y) \approx \sum_{i,j=1,...,6} u_1^{i,j } B_{i,2}(x)B_{j,2}(y) \)
\( u_2(x,y) \approx \sum_{i,j=1,...,6} u_2^{i,j } B_{i,2}(x)B_{j,2}(y) \)
oraz ciśnienia
\( p(x,y) \approx \sum_{i,j=1,...,6} p^{i,j } B_{i,2}(x)B_{j,2}(y) \)
Do testowania używamy również kombinacji liniowych funkcji B-spline
\( v_1(x,y) \approx \sum_{k,l=1,...,6} v_1^{k,l } B_{k,2}(x)B_{l,2}(y) \)
\( v_2(x,y) \approx \sum_{k,l=1,...,6} v_2^{k,l } B_{k,2}(x)B_{l,2}(y) \)
oraz ciśnienia
\( q(x,y) \approx \sum_{k,l=1,...,6} q^{k,l } B_{k,2}(x)B_{l,2}(y) \)
Przyjmujemy teraz najpierw funkcje bazowe \( ({\bf v },q) = (v_1,0,0)=(B_{k,2}(x)B_{l,2}(y),0,0) \textrm{ dla }k=l,...,6 \) oraz wstawiamy do równania . Dostajemy macierze
\( A1_{i,j=1,...,6;k,l=1,...,6} \int_{\Omega} \frac{B_{i,2}(x)}{\partial x}B_{j,2}(y)\frac{\partial B_{k,2}(x) }{\partial x}B_{l,2 }(y)dxdy + \int_{\Omega} B_{i,2}(x)\frac{\partial B_{j,2}(y)}{\partial y}B_{k,2}(x)\frac{\partial B_{l,2}(y)}{\partial y }dxdy \)
\( B1_{i,j=1,...,6;k,l=1,...,6} = \int_{\Omega} \frac{\partial B_{i,2}(x) }{\partial x } B_{j,2}(y) B_{k,2}(x)B_{l,2}(y) d,6) \)
Następnie przyjmujemy funkcje bazowe \( ({\bf v },q) = (0,v_1,0)=(0,B_{k,2}(x)B_{l,2}(y),0) \textrm{ dla }k=l,...,6 \) oraz wstawiamy do równania . Dostajemy macierze
\( A2_{i,j=1,...,6;k,l=1,...,6}=\int_{\Omega} \frac{\partial B_{i,2}(x)}{\partial x} B_{j,2}(y)\frac{\partial B_{k,2}(x)}{\partial x}B_{l,2}(y)dxdy + \int_{\Omega} B_{i,2}(x)\frac{\partial B_{j,2}(y)}{\partial y }B_{k,2}(y)\frac{\partial B_{l,2}(x)}{\partial y }dxdy \)
\( B2_{i,j=1,...,6;k,l=1,...,6} =\int_{\Omega} B_{i,2}(x)B_{j,2}(y) B_{k,2}(x)\frac{\partial B_{l,2}(y)}{\partial y } dxdy \)
Dodatkowo przyjmując funkcje bazowe \( ({\bf v },q) = (0,0,q)=(0,0,B_{k,2}(x)B_{l,2}(y)) \textrm{ dla }k=l,...,6 \) oraz wstawiamy do równania . Dostajemy macierze
\( Q1_{i,j=1,...,6;k,l=1,...,6} =\int_{\Omega } \frac{\partial B_{i,2}(x)}{\partial x} B_{j,2}(y) B_{k,2}(x) B_{l,2}(y) dxdy \)
\( Q2_{i,j=1,...,6;k,l=1,...,6} =\int_{\Omega} B_{i,2}(x)\frac{\partial B_{j,2}(y)}{\partial y}B_{k,2 }(x)B_{l,2 }(y) dxdy \)
oraz prawą stronę
\( F1_{k,l=1,...,6}= \int_{\Omega} \frac{\partial h(x,y)}{\partial x}\frac{\partial B_{k,2}(x)}{\partial x } B_{l,2}(y)dxdy +\int_{\partial \Omega} \frac{\partial h(x,y)}{\partial n } B_{k,2}(x)B_{l,2}(y)dS \)
\( F2_{k,l=1,...,6}= \int_{\Omega} \frac{\partial h(x,y)}{\partial y}B_{k,2}(x)\frac{\partial B_{l,2}(y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial h(x,y)}{\partial n } B_{k,2}(x)B_{l,2}(y)dS \)
Nasz układ równań wygląda następująco
\( \left( \begin{matrix} A1 & 0 & B1 \\ 0 & A2 & B2 \\ Q1 & Q2 & 0 \end{matrix} \right) \left( \begin{matrix} u_1 \\ u_2 \\ p \end{matrix} \right) =\left( \begin{matrix} F1 \\ F2 \\ 0 \end{matrix} \right) \)
Dokładne rozwiązanie równania Stokesa wymaga dodatkowej stabilizacji równania. W kolejnych modułach omówione zostaną dwie metody stabilizacji, jedna oparta na metodzie minimalizacji reziduum, druga oparta na metodzie DG (Discontinuous Galerkin).